Stable structures with high topological charge in nonlinear photonic quasicrystals 



K. J. H. Law,^'^ Avadh Saxena,^ P. G. Kevrekidis,^ and A. R. Bishop^ 

^ Warwick Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK 
^Department of Mathematics and Statistics, University of Massachusetts, Amherst MA 01003-4515, USA 
^Theoretical Division and Center for Nonlinear Studies, 
Los Alamos National Laboratory, Los Alamos, NM 87545, USA 

Stable vortices with topological charge of 3 and 4 are examined numerically and analytically in 
photonic quasicrystals created by interference of 5 as well as 8 beams, for cubic as well as saturable 
nonlinearities . Direct numerical simulations corroborate the analytical and numerical linear stability 
analysis predictions for such experimentally realizable structures. 
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Introduction. The study of vortices has been a princi- 
pal theme of interest in dispersive nonlinear systems with 
applications including, among others, Bose-Einstein con- 
densates (BEC), and nonlinear optical media [1-3]. More 
recently, such states have been studied in settings with 
some discrete spatial symmetry i.e., nonlinear lattices. 
There, the notion of "discrete vortices" Q arose and was 
subsequently intensely studied in both discrete and quasi- 
continuum media; see e.g. [s', ^ for relevant reviews. 
This led to the experimental realization of unit-charge 
{S = 1) coherent structures in saturably nonlinear pho- 
torefractive media [such as SBN:75(Sro.75Bao.25Nb206)] 
in 0, HI, and the exploration of higher charge {S = 2) 
ones in square and hexagonal/ honeycomb lattice settings 
[9^. A multipole soliton necklace of out-of-phase neigh- 
boring lobes in a square lattice was identified experimen- 
tally and theoretically in [10] from initial condition of a 
wide 5* = 4 gaussian beam. 

While regular lattices have been mostly studied 
more recently experimental developments have enabled 
the study of photonic quasicrystals in photorefractive me- 
dia , and have spurred a correspondingly intense theo- 
retical activity . We also note that recent experiments 
have emerged on non-square optical lattices for ultracold 
atoms in the BEC case [13]. It is then natural to expect 
that quasi-crystals are well within experimental reach in 
this regard, as well. 

Motivated by these developments, we illustrate the 
unique ability of such lattices (with saturable or cubic 
nonlinearity) to sustain stable vortices of higher topolog- 
ical charge, such as 5* = 3 and 5 = 4. Direct numerical 
simulations reveal the robustness of such modes. On the 
other hand, perhaps counter-intuitively (but as can be 
analytically predicted), lower charge vortices are found 
to be unstable, and this instability is also dynamically 
monitored. 

Theoretical Setup. We introduce the following non- 
dimensionalized evolution equation: 
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The (saturable) photorefractive nonlinearity is 
F(|/7p) = -1/(1 + |/7p) + 1, where U is the slowly 
varying amplitude of a probe beam normalized by the 
dark irradiance of the crystal Id [3, 14], and V represents 
modulation of the refractive index from interfering 
linearly propagating waves normal to the probe beam. 
In a Kerr medium the nonlinearity reads F(\U\^) = 
and this case also includes the interpretation of U as 
a mean-field wavefunction of an atomic Bose-Einstein 
condensate [15], while the potential V is either modu- 
lation of the refractive index in the former case or an 
externally applied field in the latter. 
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FIG. 1: (Color online) The stable 5* = 4 vortex in a quasi- 
crystal lattice of A/" = 5 and with a saturable nonlinearity. 
The profile and phase are depicted in panels [a(.i)], the linear 
spectrum in panel (b), Fourier spectrum in the inset panel 
(b.i), and continuation of the power, P = f \U\'^dx, as a func- 
tion of the propagation constant, /3, in panel (c). The N = 5 
lattice is depicted in (d). 



The potential V is taken to be of the form E/ (l+/(x)), 
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where /(x) 
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paradigm, this is the optical lattice intensity function 
formed by N interfering beams in the principal directions 
hj with periodicity 27r//c. We will consider the cases of 



. In the photorefractive 
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N = 5 and N = S. Here 1 is the lattice peak intensity, 
z is the propagation constant, x = (x, y) are transverse 
distances, k = 27r/5 is the wavenumber of the lattice, and 
£^ = 5 is proportional to the external voltage. Recently, 
such a setting has been explored theoretically for positive 
lattice solitons [HI, |l6| , but we extend the considerations 
here to vortex solutions. 




FIG. 2: (Color online) Panels (a-c) are the same as Fig. [T] 
except for a cubic nonlinearity. Panel (d) shows the growth 
rate, or maxA[i?e(A)]. The insets, (c.i) and (d.i) depict the 
profile and linear spectra, respectively, of the highly unstable 
solution indicated by a red square on the branches in (c,d), 
which collides with the main branch and disappears in a sad- 
dle node bifurcation close to the phonon band edge. 

The possible charge, 5, of vortices (the number of 27r 
phase shifts across a discrete contour comprising the so- 
lution) is bounded by the symmetry of the lattice [l7| . 
A lattice with n-fold symmetry has natural contours of 
2n sites. Hence, taking into account the degeneracy of 
vortex anti- vortex pairs {S^ —S}^ one has < 5 < n, 
with the cases of S = 0, n being the trivial flux cases 
of in-phase and out-of-phase neighboring lobes, respec- 
tively. The quasi-crystal with = 5 has n = 5, while 
for = 8, n = 4. Hence, the highest possible charge, 
5* = n — 1, is 5 = 4 for the case of A" = 5 and 5 = 3 for 
N = S. 

Considering the quasi-one-dimensional contour of ex- 
cited sites (depending on the respective amplitudes of 
the lattice and the probe field), and within the context 
of coupled mode theory [l8| in which the probe field is 
expanded in Wannier functions [l9|, one can obtain in- 
sights about the stability of the vortices within the frame- 
work of a discrete Nonlinear Schrodinger equation [6], 

iUn = —£{Un-\-l -\-Un-l " 2'Un) " I'^nP'^n- In that COUtCXt 

and based on either the modulational instability of jl8| , 
or through empirical numerical testing [17] or, more rig- 
orously, via Lyapunov- Schmidt perturbative expansions 
around the so-called anti- continuum (AC) limit of zero 
couphng {e = 0) [20], it is known that lobes which are 



phase-separated by greater than 7r/2 are stable next to 
each other, while those separated by less than 7r/2 are 
unstable. A simple intuitive argument for this situation 
is that the effective potential which out-of-phase neigh- 
boring nodes exert on one another through the focusing 
non-linearity is repulsive, and, hence, remain localized 
in their respective separate wells. On the other hand, 
if the neighbors are in-phase, then the effective neigh- 
boring potentials are attractive and hence the solution 
is unstable to remaining localized in separate wells. The 
possible relative phases interpolate between these cases, 
with 7r/2 being exactly in the middle. A similar discus- 
sion is used in [21] in order to justify (upon suitable phase 
variation) the existence of soliton clusters in bulk media. 
This leads to stability of the higher charged vortices for 
contours of larger numbers of nodes (see also We 
briefly review the Lyapunov- Schmidt argument. In the 
limit € ^ one can construct exact solutions of the form 
Uj = y7Ie-^-^(^^-^^)> for any arbitrary Oj G [0,27r) [2o|. 
The case we are considering is that of 6j = jSn/n. We 
linearize around the solution for £ = and the condi- 
tion for existence of solutions with e > reduces to the 
vanishing of a vector function g(^) of the phase vector 
= (^1, . . . , On)^ where 

gj = sin(^,_i - Oj) + sin(^,+i - Oj), (2) 

subject to periodic boundary conditions. This includes 
the discrete reduction of the vortex solutions for < 6* < 
n above. The fundamental contours M will have length 
\M\ = 2n, and — (/>j| = A(/) = irS/n is constant for 

ah j e M, |6>i - 0\M\ I = Ai9 and AO\M\ = mod 27r. 

For the contour M, there are \M\ eigenvalues of 
the \M\ X \M\ Jacobian Mjk = dgj/dOk of the diffeo- 
morphism given in Eq. (|2]). The eigenvalues of this 
matrix 7^ can be mapped eigenvalues of the full lin- 
earization. In particular, eigenvalues of the lineariza- 
tion, denoted Aj, are given to leading order by the re- 
lation [20] Xj = ±y^2jje. Thus, solutions are stable to 
leading order if < (so Xj G iM) and unstable if 

jj > (so A^- G M). We have jj = 4cos(A0) sin^ [0]) 
and so these cases correspond exactly to A<p > 7r/2 (or 
S > n/2) and A(/) < 7r/2 (or S < n/2). In the bound- 
ary case of A(j) = 7r/2, one needs to expand to the next 
order in the Lyapunov- Schmidt reduction. We note that 
a so-called staggering transformation along the contour, 
Uj = {—lyuj allows the above conclusions for the focus- 
ing problem to be mapped immediately to the defocus- 
ing problem (with a change in the sign of the nonlinear- 
ity). We do not consider the defocusing case further here. 
The above considerations illustrate the expectation that 
S = 3 vortices may be stable in the A^ = 5 and A^ = 8 
cases, and the 3 = 4 vortex may be stable in the A^ = 5 
case. 

Numerical Results. We now turn to numerical com- 
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put at ions. We also explore the evolution of different S 
radial gaussian beams. 
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FIG. 3: (Color online) Panels (a,b) are the same as the pre- 
vious figures for the stable S = 3 vortex in the N = 8 quasi- 
crystal lattice (c). Panel (d) is the S = 3 vortex for N = 5 
and [d.(i,ii)] are the phase and Fourier spectrum, respectively, 
of this solution. For both solutions, /3 = 3.4. 

First, we confirm the expectation of stability of the 
S = 4 vortex for saturable and cubic nonlinear it ies, over 
continuations in the semi-infinite gap (see Figs. [1] and 
[21 respectively). The profiles and phases are depicted in 
panels [a(.i)], linear spectra in panels (b), Fourier spectra 
in the inset panels (b.i), and continuations of the power, 
P = J |t/p(ix, as a function of the propagation constant, 
/3, in panels (c). The power of the solution branches dif- 
fers substantially between nonlinear it ies, and the power 
of the branch of saturable solutions approaches some res- 
onant frequency at which dP/d/3 oo and P ^ oo (see 
Fig. [11(c)). The lattice is depicted in Fig. [11(d), while 
Fig. [21 (d) shows the maximal perturbation growth rate, 
or maxA[i?e(A)], corresponding to the branches in Fig. [21 
(c). 

For the structures we consider, there is one pair of 
eigenvalues at the origin accounting for the U(l) (phase) 
invariance and the other 2n — 1 eigenvalue pairs associ- 
ated to the excited lobes all have negative energy, hence 
being candidates for instability 



22|, and are all either 



purely imaginary or purely real. If real, the instability is 
immediate, while if imaginary, instability may still arise 
due to their collision with the phonon band, resulting in 
a Hamiltonian-Hopf bifurcation and eigenvalue quartets. 
The spectral plane, with the negative energy modes indi- 
cated by red squares, for the saturable and cubic cases are 
given in panels (b) of Figs. [H and [21 respectively. Panel 
(b.ii) in Fig. [H is a closeup of the origin showing the 9 
negative energy pairs close to the origin and the one pair 
at the origin. The potential instability arising from these 
negative energy modes is prevented by their proximity 




FIG. 4: (Color online) Initial conditions (a.i,b.i) and pro- 
files at a later time (a.ii, b.ii) of the S — 4 and S — 2 radial 
Gaussian initial conditions for a saturable nonlinearity, with 
a "tight dissipation layer" (see text), D = 16. The bottom is 
similar to the top but for D — 24. 



to the origin, and distance from the phonon band. The 
expected saddle-node bifurcation [23, 24] occurs close to 
the band edge (which we computed as ~ 3.9) in which 
the main solution collides (and disappears) with an un- 
stable solution branch of a configuration with additional 
populated sites external (and in phase) to the original 
contour. 
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FIG. 5: (Color online) The dynamics of the unstable S — 2 
vortex in the case of a cubic nonlinearity. Evolution of the 
same solution with the same perturbation of random noise 
with 5% of the initial maximum amplitude of the field can 
lead to robust structures that persist for long distances (a) 
and almost immediate collapse in different trials (b). 



Next, we present results of the 5* = 3 vortex in both 
the A/" = 8 (Fig. [3l^a,b)) and AT = 5 (Fig. [31(d)) cases for 
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p = 3.4. Panel (c) depicts the N = 8 lattice, and [d.(i,ii)] 
are the phase and Fourier spectrum, respectively, of the 
solution in (d). These solutions are both stable, and 
again there is a resonance in the semi-infinite gap (not 
shown) similar to what was seen in Fig. [TJ The vortices 
for < 3 are unstable (not shown). 

To examine the potential experimental realizability 
of the above waveforms, we launch a radial Gaus- 
sian beam with topological charge 6* = 4 of the form 
^ise-{r-R) /{2h ) ^\^]^ (t^6) denoting polar coordinates, 
R = 8.5, approximately the radius of the contour, and 
6 = 1, as an initial condition into the system with sat- 
urable nonlinearity and monitor the evolution. After a 
transient period, the configuration indeed settles into an 
8 = 4: vortex contour. However, notice that this simula- 
tion has been performed with a "tight dissipation layer" 
i.e. using an extra term —iV on the right-hand side of 
Eq. [U of r = l-tanh(L)-r) with D = 16. This initially 
absorbs the shed radiation, and subsequently affects very 
little the intensity distribution. However, the phase dy- 
namics may be sensitive to the presence of such a layer: 
imposing a wider such dissipation layer with I) = 24, 
the solution actually never settles into one of constant 
charge; this topological instability effect has been ana- 
lyzed e.g. in |25l428|. Specifically, vortices may nucleate 
in the very small amplitude region and pass in and out of 
the main configuration (without affecting its intensity). 
Note that the above suggests that such effects could be 
avoided experimentally if some form of dissipation is im- 
posed. For comparison, we launch a similar initial con- 
dition with S = 2 and notice that it never settles into a 
stable configuration of fixed charge independently of the 
dissipation layer size (and despite its seemingly robust 
intensity distribution). 

Figure m (top two panels) present the initial conditions 
(a.i,b.i) and profiles for a long evolution (a.ii,b.ii) of the 
5* = 4 and S = 2 initial conditions, respectively, for sat- 
urable nonlinearity and D = 16. The charge of each fiuc- 
tuates, as power is shed and vortices nucleate in the sur- 
rounding low amplitude regions and enter and leave the 
contour as the solution traces a stationary state. How- 
ever, for the S = 4: initial condition, the field settles into 
a solution of constant charge 4 for D = 16, while for the 
S = 2 initial condition, the phase continues to fiuctu- 
ate throughout the numerical experiment. These results 
are typical in this setting. For D large (bottom panel of 
Fig. H]) the charge may never settle (topological phase 
instability). However, this does not contradict the lin- 
ear stability results (which we have confirmed separately 
for near stationary configurations). The intial condition 
^ise-{r-R) /{2b ) cos^(5^) is far from a stationary configu- 
ration and not sufficiently modulated to prevent contam- 
ination of the resulting state by radiation, although e.g. 
^io^^gi/c57r/5-(x-c.,)2-(y-c^,)2 ^-^j^ {c^j^^cyk) the ccutcr 

of one of the wells, is sufficiently localized, and with the 
latter initial condition, we observe topological stability 



(indeed without the initial turbulent fluctuating regime) 
for 5 = 4 (see Fig. |4]) but not for 5' = 2. Finally, Fig. 
[5] shows the evolution of unstable {S = 2) vortices in the 
presence of a cubic nonlinearity. The evolution depends 
sensitively on the particular initial condition. Using the 
initial condition u = U{l-\-X) with X ~ 0.05maxx[/7(x)] 
uniform over [0,1], two different particular instances can 
lead to significantly different dynamics. Either the phase 
merely reshapes as for the saturable nonlinearity, but the 
structure persists [see Fig. [5] (a)], or the solution col- 
lapses almost immediately, as seen from the maximum 
amplitude of the field in Fig. [5] (b). For larger additive 
noise, collapse seems more likely from several sample tri- 
als. The relevant mechanism involves one of the solution 
lobes exceeding a minimum collapse threshold, leading to 
an "in lobe" collapse. 

Conclusions and future directions. We have demon- 
strated numerically stable vortices of topological charge 
5* = 3 in quasi-crystals with n = 4 and 5 directions of 
symmetry and 5* = 4 with n = 5, in the cases of both 
cubic and saturable focusing nonlinearities. The nega- 
tive energy modes for these configurations remain close 
to the origin in the spectral plane, preventing collision 
with the phonon band, and can be experimentally realiz- 
able in photonic quasi-crystals in a photorefractive (or a 
Kerr) medium. This has additionally been demonstrated 
by simulation of the evolution of a radial Gaussian beam 
into such robust vortex states. This is a prime prospect 
for an immediate future experimental direction related 
to the present work. 
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